Mammal Tidal Volume

Data

Predictors: body mass interacting with habitat.

Read in data

mammal_vt <- read.csv("data/vt.csv")

Notes

log10vt is the base-10 logarithm of tidal volume in ml

Cleaning

Rename some variables and create species name from genus and species.

mammal_vt <- mammal_vt %>%
  rename(source = Source,
         order = order.corrected,
         mass.kg = body.mass..kg.) %>%
  mutate(order = str_to_title(order),
         genus = str_to_title(genus),
         animal = paste(genus, species),
         above.10kg = ifelse(mass.kg >= 10, 'Larger', 'Smaller'))

Keep only variables we will be using. And factor() “chr” variables.

mammal_vt <- mammal_vt %>%
  dplyr::select(order,
                genus,
                species, 
                common.name,
                mass.kg,
                log10.body.mass,
                log10vt,
                habitat,
                source,
                animal,
                above.10kg
         ) %>%
  mutate(across(where(is.character), factor)) %>%
  arrange(order, genus, species)

Viz

A few quick graphs just to make sure the data are looking as we expect (error checking).

my_scatter_plot <- gf_point(log10vt ~ log10.body.mass | habitat,
         color = ~order,
         # size = ~parse_number(n_animals),
         data = mammal_vt,
         alpha = 0.5) %>%
  gf_theme(legend.position = 'bottom',
           legend.title = element_text(size = 8),
           legend.text = element_text(size = 6)) %>%
  gf_theme(scale_color_viridis_d('Order')) %>%
  gf_labs(x = 'Log10(Mass (kg))',
          y = 'Log10(Vt (liter))')

my_scatter_plot

plotly::ggplotly(my_scatter_plot) %>%
  plotly::layout(legend = list(#orientation = 'h',
                               font = list(size = 6)))

GLS

Will not account for phylogeny at all in the model structure. Predictors: body mass interacting with habitat.

lm_model <- lm(log10vt ~ log10.body.mass * habitat,
                 data = mammal_vt)
tab_model(lm_model)

	log 10 vt
Predictors	Estimates	CI	p
(Intercept)	1.60	1.22 – 1.98	<0.001
log10.body.mass	0.94	0.80 – 1.08	<0.001
habitat [land]	-0.63	-1.02 – -0.25	0.002
log10.body.mass * habitat [land]	0.07	-0.07 – 0.22	0.311
Observations	50
R² / R² adjusted	0.990 / 0.990

lm_anova_results <- car::Anova(lm_model)
pander(lm_anova_results)

Anova Table (Type II tests)
	Sum Sq	Df	F value	Pr(>F)
log10.body.mass	76.68	1	2682	1.949e-42
habitat	1.376	1	48.12	1.134e-08
log10.body.mass:habitat	0.03001	1	1.05	0.3109
Residuals	1.315	46	NA	NA

Model Assessment

lm_preds <- predict(lm_model, se.fit = TRUE)

mammal_vt <- mammal_vt  %>% mutate(lm_resids = resid(lm_model),
         lm_fitted = lm_preds$fit,
         lm_fit_lo = lm_preds$fit + 1.96*lm_preds$se.fit,
         lm_fit_hi = lm_preds$fit - 1.96*lm_preds$se.fit)
gf_point(lm_resids ~ lm_fitted, data = mammal_vt)

s245::gf_acf(~lm_model)

gf_dhistogram(~lm_resids, data = mammal_vt,
              bins = 20) %>%
  gf_fitdistr()

Model Predictions

Any predictors not shown in a plot are held constant at their mean or most common value.

gf_line(lm_fitted ~ log10.body.mass,
         color = ~habitat,
         data = mammal_vt) %>%
  gf_ribbon(lm_fit_lo + lm_fit_hi ~ log10.body.mass,
            color = ~habitat, fill = ~habitat) %>% 
  gf_theme(scale_color_manual(values = my_colors)) %>%
  gf_theme(scale_fill_manual(values = my_colors))

gf_point(log10vt ~ lm_fitted, data = mammal_vt,
         alpha = 0.1) %>%
  gf_labs(y = 'Observed log10(Vt)',
          x = 'Model-predicted log10(Vt)',
          title = 'Linear Regresssion (no phylogeny)') %>%
  gf_abline(slope = 1, intercept = 0, color = 'black', linetype = 'dashed')

Mixed-effects model

Will include nested random effects of order/family/genus/species, expecting similarity of observations based on a hierarchy of phylogenetic relatedness, but not in a specified/structured way; only groupings are used, with no sense of (for example) the fact that two orders are required to be “further” apart than any other two.

re_model <- glmmTMB(log10vt ~ log10.body.mass * habitat + (1 | order/genus/species),
                 data = mammal_vt)

summary(re_model)

##  Family: gaussian  ( identity )
## Formula:          
## log10vt ~ log10.body.mass * habitat + (1 | order/genus/species)
## Data: mammal_vt
## 
##      AIC      BIC   logLik deviance df.resid 
##    -27.9    -12.6     22.0    -43.9       42 
## 
## Random effects:
## 
## Conditional model:
##  Groups                Name        Variance  Std.Dev.
##  species:(genus:order) (Intercept) 0.0108807 0.10431 
##  genus:order           (Intercept) 0.0069420 0.08332 
##  order                 (Intercept) 0.0004543 0.02132 
##  Residual                          0.0073999 0.08602 
## Number of obs: 50, groups:  
## species:(genus:order), 47; genus:order, 42; order, 9
## 
## Dispersion estimate for gaussian family (sigma^2): 0.0074 
## 
## Conditional model:
##                             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                  1.61392    0.19456   8.295  < 2e-16 ***
## log10.body.mass              0.93222    0.07429  12.548  < 2e-16 ***
## habitatland                 -0.65494    0.19922  -3.288  0.00101 ** 
## log10.body.mass:habitatland  0.07660    0.07244   1.057  0.29036    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

tab_model(re_model)

	log 10 vt
Predictors	Estimates	CI	p
(Intercept)	1.61	1.23 – 2.00	<0.001
log10.body.mass	0.93	0.79 – 1.08	<0.001
habitat [land]	-0.65	-1.05 – -0.26	0.001
log10.body.mass * habitat [land]	0.08	-0.07 – 0.22	0.290
Random Effects
σ²	0.01
τ₀₀ _{species:(genus:order)}	0.01
τ₀₀ _genus:order	0.01
τ₀₀ _order	0.00
ICC	0.71
N _species	47
N _genus	42
N _order	9
Observations	50
Marginal R² / Conditional R²	0.991 / 0.997

re_anova_results <- car::Anova(re_model)
pander(re_anova_results)

Analysis of Deviance Table (Type II Wald chisquare tests)
	Chisq	Df	Pr(>Chisq)
log10.body.mass	1594	1	0
habitat	40.93	1	1.58e-10
log10.body.mass:habitat	1.118	1	0.2904

Model Assessment

re_ave_preds <- predict(re_model, 
                    se.fit = TRUE,
                    re.form = ~0)
re_ind_preds <- predict(re_model,
                        se.fit = TRUE,
                        re.form = NULL)
mammal_vt <- mammal_vt %>%
  mutate(re_resids = resid(re_model),
         re_ind_fitted = re_ind_preds$fit,
         re_ave_fitted = re_ave_preds$fit,
         re_ave_lo = re_ave_preds$fit - 1.96*re_ave_preds$se.fit,
         re_ave_hi = re_ave_preds$fit + 1.96*re_ave_preds$se.fit)

# save fitted model and data
saveRDS(mammal_vt, 'fitted-models/vt-data.RDS')
saveRDS(re_model, 'fitted-models/vt-re-model.RDS')


gf_point(re_resids ~ re_ind_fitted, data = mammal_vt)

acf(resid(re_model))

#gf_acf(~re_model)
gf_dhistogram(~re_resids, data = mammal_vt) %>%
  gf_fitdistr()

Predictions from Model

gf_line(re_ave_fitted ~ log10.body.mass,
         color = ~habitat,
         data = mammal_vt) %>%
  gf_ribbon(re_ave_lo + re_ave_hi ~ log10.body.mass,
            color = ~habitat, fill = ~habitat) %>% 
  gf_theme(scale_color_manual(values = my_colors)) %>%
  gf_theme(scale_fill_manual(values = my_colors))

gf_point(log10vt ~ re_ind_fitted, data = mammal_vt,
         alpha = 0.1) %>%
  gf_labs(y = 'Observed log10(Vt)',
          x = 'Model-predicted, Species-specific log10(VT)',
          title = 'Mixed-effects Model (RE of Order/Genus)') %>%
  gf_abline(slope = 1, intercept = 0, color = 'black', linetype = 'dashed')

The graph above is a bit “cheating” as we have a random effect of species, but there are only 1-2 measurements for most of the species (nearly guaranteeing that our estimates will be nearly perfect). What if we also check out the predictions accounting for the modeled effects of Order and Genus, but predicting for the “average” species in each Genus?

no_species <- mammal_vt %>%
  mutate(species = NA)
re_genus_level_preds <- predict(re_model, 
                    se.fit = TRUE,
                    re.form = NULL,
                    newdata = no_species)

gf_point(log10vt ~ re_genus_level_preds$fit, data = mammal_vt,
         alpha = 0.1) %>%
  gf_labs(y = 'Observed log10(Vt)',
          x = 'Model-predicted, Genus-specific log10(Vt)',
          title = 'Mixed-effects Model (RE of Order/Genus)') %>%
  gf_abline(slope = 1, intercept = 0, color = 'black', linetype = 'dashed')

PGLS

Read in Tree Data

#Read in the trees from Upham et al

tree_path <- paste("data/upham-trees/Completed_5911sp_topoCons_FBDasZhouEtAl")
tree_files <- list.files(tree_path)

all_trees <- list()

for (i in 1:length(tree_files)){
  all_trees[[i]] <- read.tree(paste0(tree_path, '/',
                                     tree_files[i]))
  if (i == 1){
    treeset <- all_trees[[i]]
  }else{
    treeset <- c(treeset, all_trees[[i]])
  }
}

all_tip_labels <- purrr::map(treeset, "tip.label")
all_tip_labels <- purrr::map(all_tip_labels,
                             function(x) 
                               stringr::str_replace_all(x, pattern = '_',
                                                      replacement = ' '))

# get list of species that are in ALL the trees
for (t in 1:length(all_tip_labels)){
  if (t == 1){
    tip_labs <- all_tip_labels[[t]]
  }else{
    tip_labs <- intersect(tip_labs, all_tip_labels[[t]])
  }
}

                            
# keep only the species in mammal_bmr that are in all the trees
# on 4/14 this removes 1 species.
pgls_data <- mammal_vt %>%
  filter(animal %in% tip_labs) %>%
  droplevels()


taxonomy <- read_csv('data/upham-trees/taxonomy_mamPhy_5911species.csv')

Fit models, one model for every tree in our list.

pgls_models <- list()

for (t in c(1:length(treeset))){
  # make sure there is only one row of data per species (seems dubious??)
  pgls_rep_data <- pgls_data %>% 
    group_by(animal) %>%
    sample_n(1) %>%
    ungroup
  
  
  #Reduce the tree to only include those species in the data set
  refit_tree <- treeset[[t]]
  refit_tree$tip.label <- str_replace_all(refit_tree$tip.label, '_', ' ')
  refit_tree <- drop.tip(refit_tree, 
                         setdiff(refit_tree$tip.label, 
                                 levels(pgls_rep_data %>% pull(animal))))
  
  #Order the data set so that it is in the same order as the tip labels of the tree
  pgls_rep_data <- left_join(data.frame(tree.tip.label = refit_tree$tip.label),
                             pgls_rep_data,
                             by = c('tree.tip.label' = 'animal'),
                             keep = TRUE)
  
  # fit the model
 pgls_models[[t]] <- tryCatch({
    fittd <- gls(log10vt ~ log10.body.mass * habitat,
                   correlation = corPagel(value = 0.8, 
                                          phy = refit_tree, 
                                          fixed = FALSE, 
                                          form = ~animal),
                   data = pgls_rep_data)
    },
  error = function(cond){
    message(paste('PGLS fit failed for tree', t))
    return(NULL)
  }
  )
}
pgls_models <- pgls_models[!sapply(pgls_models, is.null)]

Note: we tried to fit 101 PGLS models (each with a different tree); of these, model fitting failed for 3.

Combine the many PGLS model runs together into one summary combined model according to Rubin’s rule.

# as.mira takes the list of models and create an object to be used by the mice package
pgls_mira <- as.mira(pgls_models)  
# # pool summarise the models using Rubin's rule corrected for small samples
pooled_pgls   <- pool(pgls_mira)
pooled_pgls_summ <- summary(pooled_pgls, type = 'all', conf.int = TRUE)

pander(pooled_pgls_summ %>% dplyr::select(term, estimate, std.error, `2.5 %`, `97.5 %`, lambda, fmi))

Table continues below
term	estimate	std.error	2.5 %	97.5 %
(Intercept)	1.636	0.2015	1.228	2.044
log10.body.mass	0.9389	0.07037	0.7966	1.081
habitatland	-0.6481	0.1959	-1.045	-0.2513
log10.body.mass:habitatland	0.08668	0.07333	-0.0617	0.2351

lambda	fmi
0.05837	0.1046
0.02877	0.07506
0.06361	0.1098
0.03793	0.0842

pgls_anova <- Anova(pgls_mira)
pander(pgls_anova)

Analysis of Deviance Table (Type II tests) (continued below)
	F	num df	den df	missing info	riv
log10.body.mass	2634	1	73026	0.03498	0.03625
habitat	40.54	1	31676	0.05313	0.05611
log10.body.mass:habitat	1.397	1	62093	0.03793	0.03943

	Pr(>F)
log10.body.mass	0
habitat	1.95e-10
log10.body.mass:habitat	0.2372

Alternative approach: using MuMIn to do model averaging. This will not necessarily weight all of the models/trees equally?

pgls_avg <- model.avg(pgls_models, 
                      rank = function(x) 1)

This gives a model that we can more easily make predictions from. The other way (with pool()) is better for getting the ANOVA results.

For the PGLS model, we also want to extract the estimate of \(\Lambda\), which tells us about how the phylogeny is affecting the correlation structure.

Previous approach was to take the mean of the estimates of \(\Lambda\) from all the individual fitted models.

lambda <- mean(unlist(purrr::map(pgls_models, function(x) x$modelStruct$corStruct)))
lambda

## [1] 0.07483631

According to this simple method our estimate is \(\hat{\Lambda} =\) 0.075.

Model Assessment

It’s not really clear how to even approach this, since we don’t really any longer expect the residuals to be normal or independent. And based on the data we know there’s not a huge issue with linearity. So…ok?

Mammal Tidal Volume

Fahlman & DeRuiter Labs

2021-07-27

Data

Read in data

Notes

Cleaning

Viz

GLS

Model Assessment

Model Predictions

Mixed-effects model

Model Assessment

Predictions from Model

PGLS

Read in Tree Data

Model Assessment